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Abstract 

High-temperature bivariate expansions have been derived for the two-spin correlation-function 
in a variety of classical lattice XY (planar rotator) models in which spatially isotropic interactions 
among first-neighbor spins compete with spatially isotropic or anisotropic (in particular uniaxial) 
interactions among next-to-nearest-neighbor spins. The expansions, calculated for cubic lattices of 
dimension d = 1, 2 and 3, are expressed in terms of the two variables K\ = J\/kT and Ki = J^/kT, 
where J± and J2 are the nearest-neighbor and the next-to-nearest-neighbor exchange couplings, re- 
spectively. This report deals in particular with the properties of the d = 3 uniaxial XY model 
(ANNNXY model) for which the bivariate expansions have been computed through the 18-th or- 
der, thus extending by 12 orders the results so far available and making a study of this model 
possible over a wide range of values of the competition parameter R = J%l J\. Universality with 
respect to R on the critical line separating the para- and the ferro-magnetic phases can be verified, 
and at the same time the very accurate determination 7 = 1.3177(5) and v = 0.6726(8) of the 
critical exponents of the susceptibility and of the correlation-length, in the three-dimensional XY 
universality class, can be achieved. For the exponents at the multi-critical (m,d,N) = (1,3,2) 
Lifshitz point the estimates 7/ = 1.535(25), u±_ = 0.805(15) and u» = 0.40(3) are obtained. Fi- 
nally, the susceptibility exponent is estimated along the boundary between the disordered and the 
modulated phases. 

PACS numbers: 05.50.+q, 05.70.Jk, 64.60.Kw 
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I. INTRODUCTION 



Bivariate high-temperature(HT) and, in some cases, low-temperature(LT) series expan- 
sions have been derived in the last four decades only for very few lattice spin models with 
interactions extending beyond nearest-neighbor(nn) sites. In principle, undertaking such cal- 
culations should be of permanent interest, because the analytic approximations based on the 
series coefficients can make a large part of the (at least) bidimensional interaction-parameter 
space easily accessible to analysis above (respectively below) the transition temperature. In 
practice, however, the presence of the next-to-nearest neighbor (nnn) interactions makes the 
derivation of adequately long expansions by the conventional graph techniques, very labori- 
ous. As a consequence, only data at relatively low orders and therefore of limited use, have 
been so far available for a handful of non-trivial models and the less tedious and more flex- 
ible approach to the study of these systems by stochastic simulations has largely prevailed 
in the literature, despite its forced limitation to a coarse-grained survey of the interaction- 
parameter space and the convergence problems often met with in particular regions of the 
phase diagrams. 

The earliest studies of short two-variable^ (or, in some cases, even three-variable^) series, 
particularly in the case of iV-vector spin systems, helped to substantiate and qualify the 
critical universality hypothesis in its statement concerning the independence of the critical 
exponents (and of the other universal quantities) on the range of the interaction. Decisive 
progress in the study of this property was achieved much later with the advent of simulation 
algorithms optimized^ for long-range interactions. It was however already known that, 
when the nnn interactions compete with the nn interactions, frustration occurs (in absence 
of disorder) and, as initially shown in the mean- field approximation 4,5 , produces "special 
effects" 6 . The physically interesting new features include : 

a) the formation of spatially modulated phases, i.e. spin configurations in which 
the order parameter varies periodically in space with a characteristic modulation 
wave-vector depending on the temperature and the ratio of the competing exchange 
coupling^^^^i^^'i 4 -^^^^; 

b) the occurrence of a special multi-critical point^, called Lifshitz point (LP), at which 
the HT disordered phase meets both the LT spatially-uniform ordered phase and the LT 
spatially-modulated ordered phase(s). 

Competing interactions and LP's are present in a variety of magnetic, ferroelectric, poly- 
meric, liquid crystal systems, in microemulsion models etc., which sometimes have not yet 
been studied in complete detail theoretically^^^^^^^ 1 ^^ 1 ^ or experimentally^ and 
continue to be actively explored. It is also worth to mention that within the lattice approach 
to Euclidean quantum field theory there is a continuing interest into models of the same 
or analogous structure^ 1 ^, which are expected to show a faster approach to the continuum 
limit. 

II. THE SPIN MODELS 

We begin with a general description of several systems of iV-vector spins with nn and nnn 
interactions on simple-cubic lattices of spatial dimensions d — 1,2 and 3, in zero external 
field, for which we have computed from the beginning or extended the bivariate HT series. 
Then we shall discuss a first brief analysis of a small sample of the large body of data so 
far accumulated, referring to a particular three-dimensional system of spins with N = 2 
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components, i.e. to a classical XY (or planar rotator) model. The case of general N will be 
the subject of forthcoming work. 

We have derived bivariate HT expansions for: 

i) a class of (^-dimensional models with isotropic interactions J\ among nn spins and 
isotropic (or anisotropic) interactions J2 among nnn spins separated by two lattice spacings 
along m < d lattice axes. They are described by the Hamiltonian: 

H nnn {v] = -2Jx ■ ~ 2^2 ■ v{r') (1) 

nn nnn 

where we have denoted by v(r) a iV-component classical spin vector of unit length situated 
at the lattice site r. In eq.([T]) the first sum, extended to nn spins, describes the interactions 
among the spin at f and the spins at the sites r' = r + Xi with Xi a unit lattice vector in 
the positive x\ direction. The second sum describes the interactions among the spin at r 
and the spins at the sites r ' = r + 2£i with % = 1, .., d in the isotropic d — axial case (i.e. in 
which spins separated by two spacings interact along all lattice axes), whereas i takes only 
the values i = d — m + 1, .., din the anisotropic m — axial case, with m < d. In this report 
we study in detail only the most interesting case of the m = 1 model in d = 3, which is 
sometimes denoted also as the three-dimensional ANNNXY model. 

ii) a class of models with isotropic or anisotropic interactions J\ among nn spins and 
J 2 among (geometric) second-neighbor (sn) spins (sometimes also called Ji — model or 
model with crossing bonds) described by the following Hamiltonian: 

H sn {v} = -2J x ^ti{f) ■ ') - 24j2v(r) ■ v(r') (2) 

nn sn 

The first sum in eq.([2]) has the same meaning as in eq.([T]), while the second sum extends to 
sn spins and describes the coupling of the spin at r with the spins at the sites r ' = r + Xi + xj 
and r ' = r + x,{ — xj with % < j = 1, .., d (i.e. the sn interaction acts along the diagonals of 
the elementary plaquettes). 

In eq.fll]) we have denoted by J\ and J2 the nn and the nnn exchange interaction constants 
respectively. We shall denote by R = J2/J1 their ratio, which measures the degree of 
competition of the couplings and therefore is usually called competition parameter. In eq. 
((SJ) J' 2 denotes the sn interaction constant and we set R! = J' 2 / J\. Our HT expansions are 
expressed in terms of the two variables Ki = (3J\ and K 2 = PJ2 {K' 2 = j3J 2 ), with (3 = 1/kT, 
k the Boltzmann constant and T the temperature. 

By the same techniques used in this paper, expansions can be derived also for a much 
wider variety of systems, but here we shall not be concerned with this possibility. 

In the case of the 1 — axial XY model in d = 3 studied in this report, our series extend 
through the 18-th order (altogether 171 nonzero coefficients in the case of the susceptibility) 
the existing 25 sixth-order results (21 nonzero coefficients, respectively). The extensions 
obtained for the other models described in i) and ii) are comparable. 

The complete set of our series coefficients is too extensive to print here. We shall however 
upload in the hep-lat archive a separate report containing a large sample of these data, 
including the correlation functions between nn and sn (or nnn) spins, the energy density, 
the susceptibility and a few correlation moments for the models we have studied. 

It is fair to note that, in spite of their high computational complexity and cost, our 
bivariate series still reach a length which may still be considered only "moderate" with 
respect to the current best standards for univariate HT expansions. For example, as far 
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as the accuracy of the numerical output of their analysis is concerned, they are not yet 
comparable with the 25-th order univariate series computed in the case of three-dimensional 
scalar spin systems with nn interactions^ only, (such as the Ising model with generic spin 
or the lattice Euclidean scalar field), or with the 26-th order series derived 2 ^ in the case of 
the two-dimensional XY model with nn interactions on the square lattice. However, if the 
comparison is limited to the bivariate series calculated until now for the systems i) and ii) 
(see below), our expansions seem to be already non-trivial enough to justify an update of 
the few existing HT studies or a first analysis of the data so far unavailable. 

We have studied the expansions of observables defined in terms of two-spin correlations, 
as functions of K\ at fixed values of R, using, for the moment, only the established single- 
variable methods 28 of series analysis, namely Pade approximants (PA) or inhomogeneous 
differential approximants (DA), and, adopting protocols of analysis well tested in earlier 
papers^ 2 ^ 3 ^^ 2 -, we have indeed produced reasonably accurate results which can comple- 
ment or improve those from different approaches. The choice of single variable methods of 
numerical analysis, which might perhaps be considered partly responsible of the limited accu- 
racy of our numerical results, was dictated only by simplicity, but it is likely that now the HT 
series are long enough to deserve the additional effort of an analysis by tools more powerful 
and better suited to describe bivariate critical behavior, such as the partial-differential 2 ^ 3 - 
approximants or even some simpler two- variable generalization of PAs 34 . 

For N > 1, the lower critical dimension of the so called (m, d, A) Lifshitz point occurring 
in a m — axial A- vector system in d dimensions, is di(m) = 2 + m/2, while the upper critical 
dimensional^ 3 - is d u (m) = 4 + m/2. Therefore one has to expect that only in dimension 
d > 2, for the 1 — axial model, a LP will show up at T^p > with a non-classical critical 
behavior. Thus we can use our HT series for the three-dimensional ANNNXY model to 
locate its LP and to obtain estimates of its critical exponents, for which some predictions 
from other approximation methods already exist. 

More generally, it is also of interest to analyze the behavior of the series as functions of R 
all along the critical line separating the paramagnetic from the ordered phases. In particu- 
lar, along the critical line between the disordered and the ferromagnetic phases we can take 
advantage in our analysis of the old suggestions^ that the accuracy in the determination 
of universal critical parameters of a given A-vector model, such as critical exponents and 
universal amplitude ratios, can be significantly improved by extending the analysis to a one- 
parameter family of models belonging to the same universality class. One should then simply 
"tune" the family parameter in order to minimize or, if possible, to suppress the amplitudes 
of the (non-universal) leading non-analytic corrections to scaling in the observables used 
to evaluate the critical parameters. This procedure can be implemented by a Monte Carlo 
(MC) method onlyS, by HT -series assisted MC 38 , or by HT series only 2 ^, and has been 
successfully applied in a variety of cases involving nn and local spin interactions, so far with 
a single exceptional for the spin- 1/2 Ising model on the simple-cubic lattice. In this case, a 
high-precision MC analysis showed that turning on an isotropic ferromagnetic coupling be- 
yond nn results into a considerable decrease of the leading correction-to-scaling amplitudes. 
Unfortunately, this simulation study had to be restricted to a single well-guessed value of the 
nnn coupling, because "tuning" the additional (irrelevant) interaction to search for its best 
value, was too time-consuming. Now, in the case of the A— vector model with A = 2, our 
expansions are sufficiently extended to produce accurate analytic approximations enabling 
us to determine in a straightforward way the optimal value of the irrelevant nnn coupling 
within a fair approximation. 
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The LP can be located accurately and the exponents of the susceptibility and of the 
transverse correlation-length can be determined with fair precision, while the exponent of 
the parallel correlation-length can be somewhat less accurately estimated. We have also 
examined the critical behavior of the susceptibility along the branch of the critical line 
between the disordered and the modulated phases. In this case the results are significantly 
less accurate and somewhat puzzling: in particular it is not clear whether the transition is 
(weakly) first-order or second-order and, if this is the case, to which universality class it 
belongs. 

Finally, we can point out another valuable use of our expansions as a guide for possibly 
more detailed MC investigations of models in the classes i) and ii), which might focus on 
specific points of the parameter space. Even more simply, our expansions can also serve as 
a realistic test-ground for techniques of analysis and re-summation of multivariate series. If 
nothing else, our calculations provide, for a variety of models, an initial set of HT reference 
data large enough to be a stringent constraint in the validation of future series extensions. 

For completeness and in order to put our work into perspective, it is convenient to list 
the main existing HT results for the systems i) and ii) and to mention a few studies related 
to the subject of this paper, but using different techniques. 

The earliest bivariate HT series investigation 8 - of the LP in the iV-vector model with 
nn and uniaxial nnn interactions, was mainly devoted to the N = 1 case (namely the 
spin-1/2 Ising model). In three dimensions, series of order 8,6 and 5 were derived for 
N — 1,2 and 3, respectively. However, only the Ising HT series were considered^ by 
the authors long enough to yield sufficiently reliable numerical estimates of the location 
and the exponents of the LP. Both in two and three dimensions, the N — 1 series were 
subsequently extendecP^ 1 ^ through order 11 for the J± — Jg interaction, for the 1-axial 
and for the d-axial interactions. Soon later, series for the susceptibility through the twelfth 
order were computed and analyzed^ 2 - for uniaxial models on the simple-cubic and the face- 
centered-cubic lattices. More recently, in two dimensions for the square-lattice Ising model, 
the susceptibility expansion was pushed— through the thirteenth order for the Ji — J2 and 
the 2-axial interactions. These have been until now the longest bivariate series derived for 
systems with nn and nnn interactions. 

For the N > 2 models, on the other hand, there has been no further progress in the 
HT calculations for the last three decades since the early sixth-order results 25 in three 
dimensions. The only exception to this lack of activity was a study of the 2-axial N- vector 
model (with isotropic nn and nnn interactions on the square lattice), in which the bivariate 
HT series expansion coefficients were computed and tabulated through the fifth-order—, for 
general values of N. Interest into this model, which constitutes the "Symanzik improved"— 
square- lattice formulation of the 0(7V)-symmetric non-linear cr-model in two-dimensional 
Euclidean field theory, came however from quantum field theorists. In particular, these HT 
expansions were derived as a means to infer the weak-coupling (i.e. LT) properties of the 
"improved" non-linear cx-model for N > 3. 

A variety of MC investigations of models with nnn interactions can also be found in the 
literature: most of them are concerned with the N = 1 caiS Q§.i2£.£ldl.d2.&kM.dk j and some, more 
directly related to our work, with the N = 2 case. In particular, the uniaxial XY model 
in three dimensions was simulated in Refs. |47I |. its phase diagram was sketched out and 
the location of the LP was estimated along with the exponents of the susceptibility and the 
magnetization. Later investigations^ 1 ^ were devoted also to the J\ — J2 XY model in two 
dimensions. These studies were generally limited to only a few points in the interaction- 



5 



parameter space. 

Before a systematic renormalization group approach could be extended to cover also the 
class of systems studied here, non-trivial technical difficulties, due to the anisotropy of the 
scale invariance in the critical behavior at LP's, had to be solved. In particular, while the 
lowest order computation of all critical exponents by the e-expansion around the LP upper 
critical dimension d u {m) (with e = d u (m) — d), as well as some 0(e 2 ) results, go back to three 
decades ago&, the extension through second order has been completed^^i only recently. It 
is encouraging that the 0(e 2 ) corrections are small and decrease with increasing N, so that 
a simple truncation to order e 2 of the exponent expansions might already lead to reasonable 
approximations. Of course, further work is still needed before an accuracy comparable to 
that established for the usual critical points can be attained. 

Results consistent with the 0(e 2 ) calculations for the exponents of the Ising 1 — axial 
model in three dimensions have also been obtained^ 2 - from a truncation of the exact renor- 
malization group equation. Renormalization group discussions of the ANNNXY model have 



There has been recent progress^ also in the calculation of the leading non-trivial correc- 
tions to the large N limits for some exponents. 

The layout of the paper is the following: In Sect. II the quantities for which we have 
computed HT expansions are defined in detail, and related to the main critical parameters 
in particular to the exponents of the uniaxial LP. The Sect. Ill is devoted to a discussion of 
the numerical analysis of the series. In the Appendix, we briefly describe the non-graphical 
algorithm used to compute the HT expansions and the checks passed by our series data. 

III. THE HIGH-TEMPERATURE EXPANSIONS 

In zero field and for any dimension d, for all models of the classes i) and ii) defined on 
bipartite lattices, the free energy is an even function of J\ and therefore we can restrict our 
analyses to the case J\ > (ferromagnetic nn interaction). 

In general, for d > 1 these models exhibit three main phases: a HT paramagnetic (P) 
phase, a LT uniformly-ordered ferromagnetic (F) phase and a family of LT ordered modu- 
lated (M) phases. For d > 1, a non-trivial transition line K\ C {R) (K\ C (R')) separates the 
HT P-phase from the LT phases, whereas for d — 1, one has simply K\ C (R) = oo. For 
d > di(m), the critical line K\ C {R) is divided into two branches by a LP, a triple point 
located at a non-zero value K\ c (Rlp) and separating the P-F-transition line from the P-M 
line. 

For T = 0, the spin ordering can be determined simply by minimizing the energy. The 
ferromagnetic ground state is energetically favored over the modulated phase(s) only for R 
greater than some critical value. In the three-dimensional uniaxial case under study, the 
ground state is made of ferromagnetic layers orthogonal to the z-axis, with the relative 
orientation of the successive layers determined by the value of R. 

It is clear that, except in the d = 1 case, in which the LT region is simply shrunk to the 
border of the paramagnetic region, the zero-field HT expansions are unsuited to yield much 
more than hints on the LT structure of the phase diagram and therefore, for the purpose of 
investigating the LT phases, they have to be replaced by other methods. 

For all models of the classes i) and ii), we have derived the bivariate expansion of the 
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spin-spin correlation-function, 

C(0, x; K h R) =< v{0) ■ v(x) >, (3) 

for all values of x for which non-vanishing coefficients exist within the maximum order of 
expansion. 

The appropriate quantities to be studied in order to locate the P-F branch of the critical 
line are the ordinary susceptibility 

X (K 1 ,R) = l + J2<v(0)-v(x) > (4) 

x^O 

and the Z-th order spherical moments of the correlation function 

m®(K x , R) = < ^(0) ■ v(x) > ■ (5) 

X 

In terms of va^ 2 '(Ki, R) and xC^ij we can construct the correlation length 

e(Kt,R) = v^\K x ,R)/2d X {K u R) (6) 

When studying am — axial model with m < d and therefore with anisotropic interactions, 
it is convenient to break the d— dimensional lattice vectors x as x = (xn,x±), where x» and 
x± denote the m-dimensional and the (d — m)-dimensional components of x, respectively 
parallel and perpendicular to the directions of the nnn interaction. 

In order to study the properties of the LP and to locate the P-M transition from the 
paramagnetic to the LT modulated phase it is necessary to study also the structure function 
with respect to qu 

X (q\i ;K 1 ,R) = 1 + J2 < v(0) ■ v(x) > (7) 

x^O 

and to q± 

X (q±; K U R) = 1 + J2 < ^(0) ■ v(x) > . (8) 

x^O 

In a vicinity of the LP, the scaling behavior of the two-spin correlation function becomes 
strongly anisotropic. In particular, the asymptotic behavior for large separation of spins 
joined by a vector whose components lie entirely in the m— dimensional x\\ subspace, differs 
from the behavior in the case in which the spins are joined by a vector in the (d — Tri- 
dimensional x± subspace. It is then necessary to replace^ each one of the usual correlation 
exponents rj and v by a pair of exponents associated to the two subspaces. In the first 
case, the 'parallel correlation-exponent is usually denoted by 77/4 (or less frequently, but more 
suggestively by 7711), in the second case the transverse correlation-exponent is denoted by 
rji2 (or by 1]±). Correspondingly, in different directions two distinct correlation lengths £y 
and are observed which diverge with different exponents uu and u±, respectively. These 
exponents are related to the susceptibility exponent 7^ at the LP, by the "anisotropic scaling 
laws" 

7i = (2-7?±K = (4-77||)i/||. (9) 
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Moreover a generalized hyperscaling law is expected to hold: 



2 — ai = mv\\ + (d — m)h , ±. (10) 

The other scaling relations: a\ + 2/5; + 7^ = 2 and 7; = /3i(5i — 1) remain unchanged^. We can 
thus conclude that three independent exponents are requested to characterize the uniaxial 
LP. 

When d > d u (m) = 4 + m/2, and m = 1, the critical exponents assume the 
iV— independent mean-field values reported in Table [H 

In order to estimate the additional critical exponents characterizing a LP, we have also 
computed the l-th order "parallel moments" of the correlation function 

mf{K!,R) = ^2\x\\\ l <v{0)-v{x) > (11) 

X 

and the l-th order " perpendicular moments" 

m^(K u R) = J2 \x±\ l < 0) ■ v(x) > . (12) 

X 

Near the LP, the correlation length Q[Ki, R) in a direction within the m-dimensional 
subspace of the nnn interaction is then expressed in terms of these quantities by 

R LP ) = mjj 2) (tf ls R LP )/2m X (0; K u R LP ) ~ r(R LP y 2 ^ (13) 

and analogously, the correlation length ^\{K\,R) in a direction orthogonal to the nnn 
interaction by 

R LP ) = mf {K x , R LP )/2(d - m)x(0; K u R LP ) ~ r(R LP )- 2 ^ (14) 

with r(i?) = T/T C (R) — 1. Defining, in analogy with t(R lp ), the reduced competition ratio 
Pl P = R/Rl P — 1, a crossover exponent can be introduced to characterize the behavior of 
the reduced critical temperature t lp = T C (R)/T C (R LP ) — 1 as the LP is approached along 
the critical line: r LP ~ \plp\ 1 ^ Beside the above exponents, a related^ one: (3 q = v\\/(f) is 
associated to the LP. It describes the behavior of the magnitude of the modulation vector q 
as the LP is approached along the branch of the critical line separating the disordered and 
the modulated ordered phases: \q 2 \ ~ p^p. 

IV. THE THREE-DIMENSIONAL UNIAXIAL (ANNNXY) MODEL 
A. Universality along the P-F branch of the critical line 

For R in the range —1/4 < R < 00, the ground state of the system is ferromagnetic 
and the uniform ferromagnetic ordering persists at T > up to some inverse temperature 
Ki c (R), at which a second-order phase transition, expected to belong to the universality class 
of the three-dimensional XY model, occurs between the LT phase and the HT paramagnetic 
phase. 

In this subsection we shall mainly discuss the numerical evidence obtained from the 
analysis of our HT series that, in the ferromagnetic range of R, this transition actually 



belongs to the XY universality class. We shall then argue that, if this is the case, the 
parameter R can be exploited to determine very accurate values of the critical exponents 
for the XY universality class. 

For each value of R, we can locate the transition by analyzing the HT expansion of the 
ordinary susceptibility x(J£i, R)> whose coefficients show generally a smooth dependence on 
the order of expansion and a fast approach to their asymptotic forms, and thus are well 
suited to numerical study. 

The critical behavior of the susceptibility as r(R) — > + is expected to be 

X (K lt R) = A+{R)t{R)-^ (l + a+(R)r(RyM + 0(r(R))) (15) 

where A^(R) is the critical amplitude of the susceptibility, and a^(R) is the leading 
correction-to-scaling amplitude. We have denoted j(R) and u(R) as a priori independent, 
although we shall finally argue that, as expected, they are universal with respect to R, 
namely independent of R. 

The critical behavior of the second-moment correlation length may analogously be char- 
acterized as 

e^R) = 4(fl)r(fi)- 2 "( Ji )(l + 4(fi)r(fl) 1 '( fl ) + 0(r(]J))). (16) 

We shall re-sum the susceptibility series by inhomogeneous second-order DAs in the 
variable K\ at fixed R. In this and in the analyses that follow, we have used a set of quasi- 
diagonal DAs chosen as the approximants [k, I, m; n] with 14 < k + I + m + n < 16, namely 
those using not less than 17 series coefficients. We have taken \k — l\,\l — m\, \k — m\ < 3 
with k,l,m > 3, while 1 < n < 4. We have however always made sure that our numerical 
estimates depend only weakly on this choice. The results for the critical line obtained in 
this way (for R LP < R < 2.) are reported in Fig. [TJ Notice that we have preferred to 
plot vs R the quantity T c (R)/2 = 1 /2Ki c (R) rather than Ki c (R) itself, in order to make 
our figure immediately comparable with the figure, covering a smaller range of R, which 
appears in the MC study of Ref . (46| . On the scale of Fig. [T], the three data points (R = 
Q.,T c {R)/2 = 2.17(2)), (R = -0.25, T c {R)/2 = 1.83(2)) and (R = -0.26, T c (R)/2 = 1.82(2)) 
obtained^ in the MC study, are hardly distinguishable from our curve. The corresponding 
values determined by our series are: (R = 0.,T c (R)/2 = 2.2017(2)), (R = -0.25, T c (R)/2 = 
1.830(1)) and (R = -0.26, T c {R)/2 = 1.809(1)). Our result at R = 0. compares well with the 
estimate (R = 0.,T c (R)/2 = 2.20172(15)), obtained 58 by 21th-order HT expansions. Notice 
also that the spreads SKi c (R) of our DA estimates of the critical temperatures, which are 
smaller than those of the MC study by one order of magnitude, are invisible on the scale of 
the figure. For convenience, we have listed in Table HT1 a few numerical values of K\ C (R). 

It should be stressed that it is generally difficult to assess very accurately the real uncer- 
tainties of the results in this kind of analysis mainly because, due to the finite (and in our 
case still moderate) length of the series, the sequences of DA estimates may retain residual 
trends which call for further extrapolation, particularly so for R > 1. Whenever possible, one 
should then try to infer the size of the uncertainties also from a comparison with the results 
of approximation procedures alternative to the direct DA calculation, and thus presumably 
having different convergence rates and different mechanisms of error build-up. In our case, 
we shall eventually argue that the spread of the DA estimates of the critical temperature, 
already at these orders of HT expansion, are reasonable approximations of the uncertainties, 
at least for positive and not too large R. 
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Let us then consider an example of an alternative approach to the determination of 
Ki c (R). If the critical singularity is the nearest singularity, we can determine it also by eval- 
uating the limit of the sequence of estimators (iT lc (i2)J of K\ C (R) defined by the modified- 
ratio prescription^^: 

(K lc (R)) n = (^^) 1/4 exp[ S " + ^- 2 ] (17) 

where 

a n =(H-^)- 1 +H-^-r X )l^ (18) 

and c n (R) are the HT expansion coefficients of the susceptibility. 

This prescription has the important advantage of bringing information not only about 
Ki c (R), but also on the leading correction-to-scaling amplitude a+(R), defined by eq. ffToT) . 
a quantity which in general rules the convergence properties of any approximation method 
in the critical region. Indeed, we have observed^ that the modified-ratio sequence has the 
simple asymptotic behavior for large order n 

(K lc (R)) n = K lc {R) (l - C{1 ^}$ {R) + 0(l/n 2 )) (19) 

where C(j, u) is some known positive function of 7 and lu. 

In order to use effectively the modified-ratio method, we shall now assume that the 
exponent ui = u(R) of the leading correction to scaling is independent of R, and takes 
the valued to ~ 0.52 (an assumption which was not necessary to make in the DA method 
discussion). Then by fitting the sequence {Ki c {R)') to the simple form b\ (R) — b 2 (R) / n 1+UJ , 
we can estimate the amplitude a^(R) from the value of the parameter b 2 (R). Although 
rather long and smooth series are usually necessary^ to obtain accurate estimates by this 
method, we can observe a complete consistency between the estimates of K\ C {R) from the 
DAs and the values of bi (R) obtained by the fit (within a small multiple of their spreads, for 
positive and not too large R), so that the results cannot be distinguished from those reported 
in FigJTJ This fact also gives support to our simple fit procedure for determining b 2 (R) and, 
at the same time, it suggests that the spread of the DA estimates of the critical temperatures 
might be taken as a sound measure of their real uncertainties. In order to illustrate our fit 
procedure for determining b 2 (R), in Figj2] we have plotted vs x = l/n 1+UJ the sequence of 
modified ratios (Ki c (R)) n normalized to their extrapolated values bi(R), for a few fixed 
values of R, chosen in a vicinity of Rm- Let us add that only for R > the sequences 
(Ki c (R)) n are sufficiently smooth that these rough, but sufficiently reliable, estimates of 
b 2 (R) are feasible, while unfortunately for R < the modified-ratio sequences develop 
strong oscillations and a straightforward two-parameter fit cannot work. The important 
observation is now that the function b 2 (R) vanishes at R = Rm — 0.28(3) and therefore its 
absolute value is minimum at this point. From our estimates of b 2 (R), we can infer the sign 
and size of the deviations from the exact values which should be expected for the central DA 
estimates of Ki c (R) and of the exponents, as R varies. More precisely, we have to expect^ 
that, in the range R > Rm, where b 2 (R) and therefore a^(R) are found to be positive, 
the critical inverse temperatures and the critical exponents will be underestimated by our 
analyses, while the opposite will be observed for R < R M - This is clear from eq{19], as far as 
the critical inverse temperatures are concerned. In order to reach the same conclusion for the 
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exponents, one may either use an asymptotic formula 26 analogous to eq{T9j or consider that 
in approximate calculations some effective exponent — is generally evaluated, for example: 
7e//(^) = ~ diog(r j — 7 — uja x^ UJ w ith t small, but nonzero. In order to relate these remarks 
to the behavior of 62 (R), the absolute value of 62 {R), obtained from a fit of the four highest- 
order estimators (Ki c (R)) n in the modified-ratio sequence, is very schematically plotted vs 
R, together with the exponent estimates, in Figure [3] and in some of the following figures. 
All previous considerations also apply to the study of the HT expansion of £ 2 (Ki, R). 

In Figure [3j the critical exponents 7(-R) of the susceptibility and v±{R) of the transverse 
correlation- length (notice that ^j_(-R) and v\\(R) coincide with v(R) for R > 0) are plotted 
as functions of R along the P-F branch of the critical line. They have been computed 
both using second-order DAs biased with the critical singularity K lc (R) and, alternatively, 
also by the method of "critical-point renormalization"— . The latter method consists in 
analyzing the term-by-term divided series W(x, R) = ^2 s d s (R)/c s (R)x s , where d s (R) are 
the expansion coefficients of x*(Ki, R) (or of £ 4 (ifi, R)/Kf) and c s (R) are the coefficients of 
x{K\i R) (or of £ 2 (Ki, R)/Ki). It can be shown that, if the nearest singularity is the critical 
point, then W(x,R) is singular at x — 1 with an exponent — (l + j(R)) (resp. — (1 + 2z/(i?))), 
which can be estimated by forming DAs of W(x, R) biased to be singular at x — 1. This 
alternative determination of the exponents is of particular interest in the ranges of values 
of R in which the accuracy of the available estimates of Ki c (R) might be insufficient to 
obtain good temperature-biased estimates or where the convergence of the DAs is slow. 
It is not surprising that in a neighborhood of Rm (where the amplitude of the leading 
correction to scaling is vanishing), this method and the direct analysis of the susceptibility 
(or of the correlation length) by temperature-biased DAs yield essentially the same exponent 
estimates, while elsewhere they show some small difference. In Fig. [3j the results of both 
approximations are plotted vs R and compared with a recent high accuracy determination^ 8 - 
7 = 1.3178(2) and v = 0.67155(27) of the exponents 7 and v for the XY universality class. 
Our two approximations for 'y(R) and v{R) show a similar behavior: for R > Rm both 
lead to estimates slightly smaller than the data chosen for comparison, while the opposite 
happens for R < Rm- This is precisely what can be anticipated from our determination of 
the leading correction-to-scaling amplitude. Overall, as shown in Figure [31 for 0.1 < R < 1., 
the central values of our estimates of the critical exponents deviate from the estimates at 
R = Rm by less than 0.5%, while in the wider range 0. < R < 1.5 the deviations do not 
exceed 1%, thus indicating that our approximate results have a very weak dependence on R, 
to within a fair accuracy. Even more accurate results are obtained computing, for example 

by simplified 61 DAs, (biased with Ki c (R) and with the correction-to-scaling exponent uj), 

(2) 

the ratio of the log-derivatives of the quantities m ± (Ki, R)/ K\ and X\^U R)i which yields 
the ratio h>±(R)/ , y(R). As it is also shown in Figj3j the estimates so obtained for this ratio 
appear to be independent of R to within 0.1% along the whole interval 0. < R < 1.5, in 
which they remain quite near to the ratio of the data 3 - 8 - chosen for comparison. Of course, 
this particularly favorable result is simply due to the fact that the relative deviations of 
■y(R) and u±(R), with respect to their values at Rm, keep the same sign and a similar size 
as R varies. 

A blown-up view of part of these results is presented in FigJU where we have plotted 
our temperature-biased DA estimates of the exponents j(R), v±(R) and of their ratio, after 
normalizing them to the corresponding comparison 38 values. 

It is amusing to remark that, for positive and not too large values of R, our rough 
determination of the behavior of a x (R) also suggests a simple prescription to improve the 
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estimates of the exponents by temperature-biased DAs. We have simply to correct for 
the expected errors in the bias values of Ki c (R), computing the critical exponents by DAs 
biased with K lc (R) + 8K lc (R)/2 when R > R M + 0.03, or with K lc (R) - 8K lc {R)/2 when 
R < Rm — 0.03. We have shown in FigJ5]that this quite naive adjustment of the standard 
biasing procedure to account for the sign and size of a x (R) improves visibly the universality 
of the exponents with respect to R. At the same time, this result gives further support to our 
conjecture that the spread 8K\ C (R) of the DA estimates of Ki c (R) is a sound approximation 
of their uncertainty, provided that R is positive and not too large. 

We have also studied other indicators of universality, such as, for example, the values at 
the critical point K\ C (R) of the correlation- moment ratios 

n( m ^(K^R^jK^R) 

Q ^ q > r ' S; R) = m(r)( Ku R)^(K u R) (20) 

with p + q = r + s. As expected, they show an approximate independence on R for R > 0. 
This is shown in FigJHl where we have plotted vs R a few ratios Q(p, q; r, s; R) evaluated at 
the critical point and normalized to their value Q(p, q; r, s; Rm) at Rm- Our estimates refer 
to the cases in which (q = p = 1/2; r = 1/4, s — 3/4), or (q — 1/2, p — 1/4; r = 0, s — 3/4) 
or (q = p = 1/2; r = 0, s = 1). 

As long as R is positive and not too large, we can conclude that all these results con- 
sistently and rather convincingly indicate that, along the P-F branch of the critical line, 
the small violations of the exponent universality with respect to R, shown by our numerical 
computations, are only apparent and can be entirely ascribed to the slow convergence of 
approximation procedures still unable, at the present orders of expansion, to account fully 
for the presence of corrections to scaling. 

These evidences of universality with respect to R justify the technique^ 1 ^ that we can 
adopt in order to improve the accuracy in the determination of the exponents for the XY 
universality class. We can observe that, as R varies in the ferromagnetic range, we have an 
-R-dependent family of models all of which can be assumed to belong to the XY universality 
class, so that they share the same critical exponents, while they have generally different 
-R-dependent (namely non-universal) amplitudes of the corrections to scaling. Therefore, we 
expect that the best approximations for the universal quantities will be achieved from the 
study of the model with R = R M , because a+(i?A/) vanishes. For this particular model in 
the family also the other leading correction amplitudes of interest, for example a^(-R), must 
vanish at Rm, since the correction-amplitude ratios such as a^(R)/a^ 2 (R) are universal. 

These arguments support our belief that our exponent estimates j(Rm) = 1.3177(5), 
v(Rm) = 0.6726(8) and v(Rm) /j(Rm) = 0.5100(1) should be rated as the best possible 
determinations of the susceptibility and correlation-length exponents in the XY universality 
class, that one can extract, at the present expansion order, from our .R-dependent family of 
HT series. 

The high-precision estimates^ of the XY universality class exponents that we have com- 
pared to our results in FigJ31 were also obtained using a similar improvement procedure 
in the case of a different family of nn-interaction models for which series of order 23 are 
known. Our best estimate of the susceptibility exponent agrees well with the corresponding 
result cited for comparison, although the uncertainty of our result is sizably larger, due to 
the still moderate length of our bivariate series. On the other hand, the central value of our 
best estimate for the exponent v is somewhat larger than (while the estimate v = 0.6720(4) 
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obtained from those of 7 and of the ratio 1//7 is much closer to) the corresponding com- 
parison value. Our direct estimate of v is closer to the older MC result 62 v = 0.6723(3) [8] 
and to the result v = 0.6717(3) of the simulation of Refs.[63.64|. It is, however, much larger 
than the value v = 0.6709(1) obtained, by using the hyperscaling relation 2 — a = du, from 
the recent high-accuracy measure a = 0.0127(3) of the 4 He specific heat in a micro-gravity 
experiment 6 ^. 

In Ref. 63] recent determinations of v with increasing accuracy have been summarized into 
a useful diagram showing that the central estimates from the most recent MC simulations 
and HT series analyses are, in general, significantly larger than those obtained both from 
the renormalization group and from the cited experimental measure. This is an interesting 
remark which calls, at least, for a more accurate reassessment of the uncertainties of the 
results in the literature. As far as our HT approach is concerned, we can reasonably expect 
that an extension of the bivariate expansions by only a couple of orders would significantly 
reduce the uncertainties of our best estimates of the exponents, particularly so for the direct 
estimate of v. Presently, however, it might be less difficult to give further support to our 
arguments and make them sharper by using our series results as a guide for a high-precision 
MC simulation of the ANNNXY model at R = R M - 

Finally we must add that, unfortunately, our series for the specific heat, which is a 
very weakly singular quantity, seems to be still insufficiently long to yield an evaluation of 
comparable accuracy for the exponent a(R). We can only infer that on the P-Fline, a(R) ~ 
—0.01(2), as expected. This result is completely compatible with the bounds —0.0202 < a < 
—0.0124 obtained introducing into the hyperscaling relation the extremal values of the range 
of estimates of v reported in the recent literature. 

We can conclude our study of the critical behavior along the P-F line remarking that 
from a study of T C (R) in a vicinity of T^p, we can estimate <p — 1-00(4), since it appears 
that the curve T c = T C (R) changes from concave to convex in a small vicinity of T(Rlp). 



B. The Lifshitz point and the P-M branch of the critical line 

In the mean-field approximation Rlp = —1/4, but our HT calculation shifts this value by 
~ 10% to the modulated side of the phase diagram. For Rlp < R < 0, the critical exponents 
crossover— from the value of the XY universality class to the LP critical behavior. In this 
range of R, the amplitudes a+(R) and a^ 2 (R) appearing in eqs. f|T5|) and ffTBl) lose their leading 
role, because also higher order correction amplitudes become important, the convergence of 
our series slows down and they appear inadequate to exhibit the universality with respect 
to R of the exponents and their expected discontinuous change to the values of the LP 
universality class at R = Rlp- Thus, of course, in the crossover region, also the spreads of 
our estimates will grossly underestimate the real errors. 

Let us now recall that the eq. m (Ki,R) = implicitly defines in the R — T plane 
the so-called disorder line^, which, within the paramagnetic phase, divides a region with 
monotonically (exponentially) decaying correlations along the direction parallel to the nnn 
interaction from a region with oscillating, but still exponentially damped correlations. The 
equation defining the disorder line can be solved iteratively with respect to R to form the 
single- variable series Rdi S = R<ns{T), which finally is re-summed by DAs. The plot of the 
disorder line obtained from this series is drawn as an almost vertical dashed line in Figure 
[TJ It may be of interest to show how accurately the spin-spin correlations can be computed 
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at HT and therefore we have displayed in FigJT] the qualitative difference in their behavior 
as functions of the distance of the spins along the z-direction, on the two sides of the 
disorder line. Using our knowledge of the spin-spin correlations, we can also show that 
the LT modulated order already begins to build up in the nearly critical HT phase. This is 
suggested by FiglSJ where we have plotted vs R the values of the energy, the nn and the nnn 
spin correlations along the z-axis, calculated just above the boundary of the paramagnetic 
phase, precisely at T = 1.1T C (R). We can observe that the energy reaches a maximum near 
Rip-, where the disorder is higher and that the nn spins are positively correlated on the 
whole range of R (albeit not too strongly since T is high). On the other hand the nnn 
spins tend to be more correlated than the nn spins for R > 1.2, while as R is lowered, this 
correlation decays to become negative when R < Rlp- 

In order to locate the LP on the boundary of the paramagnetic phase, we have to recall^^ 
that for small q z 

mf(K u R) A m^(K u R) 2 m« 4, ( A',, R) > 

x( o, Ki, *)/xfe, k u R) = i + + ft-fc^ - afei^) + - (21) 

showing that the minimum at q z = 0, which characterizes x(0, Ki, R)/x{ ( lzi Ki, R) at fixed 
K\, when R is in the ferromagnetic range, changes to a local maximum as R — > Rdi S {K~i) 

(2) 

where the second-order parallel moment m (Ki, R) = 0. Thus the LP is found at the inter- 
section of the P-F branch of the critical line with the disorder line. Following this procedure 
we are led to the estimate {Rlp = —0.2733(6), T(Rlp) = 1.778(2)) of the intersection point 
between the critical locus and the disorder line. The value of Rlp obtained in this way 
is consistent with that obtained minimizing ^(0, Ki, R)/x{Qzi K\-> R) with respect to q\ as 
R — Rlp — > CT and K\ — > K\{R). Indeed, for small q 2 , we obtain from eq. (j2"Tl) that the 
position of the minimum q 2 z vs R is given by 



6m ( p (K 1 ,R) X (K 1 ,R) 
mf, 4 ^!, R)x(K 1 , R) - Gmf^Ki, Rf 



t « IaSTZ J 77T^ - r^T- 7^ ( 22 ) 



evaluated at the critical point Ki c = K lc (R) for R < Rlp- By re-summing the series 
expansion of q 2 z1 we can determine Rlp, also as the value of R at which q 2 z = 0. 

For R < 0, in the crossover region, FigsJH] and M show a steep rise of our exponent 
estimates near the LP and past it, followed by a somewhat slower decay extending through 
R ~ —0.8. A completely similar behavior of the exponent 'y(R) nearby the LP was noted in 
a HT study22 of the 3d ANNNI model. In this study, we should then assume that the values 
of the exponents at the LP are not strongly affected by the crossover and also be aware that 
they are rather sensitive to the location of the LP. Using the critical-point renormalization 
method, we can estimate 7/ = 1.55(1), while the biased DAs suggest 72 = 1.52(1), as shown in 
Fig OJ We shall take a weighted average of these values as our final estimate of the exponent: 
7; = 1.535±0.025±0.2|i?£p+0.2733|, including explicitly in our error a contribution from the 
uncertainty of Rlp- This result has a smaller uncertainty than but is completely compatible 
with the estimate 7; = 1.5(1) of the MC simulation of Ref . (47| . Notice that our estimates of 
R LP and T(R LP )/2. differ by ~ 5% from the values R LP = -0.263(2) and T LP (R)/2. ~ 1.82, 
determined by the old (sixth-order) series^ and used as an input in the MG 4 ^ 1 ^ study. Our 
estimate of the exponent is also not far from the value 7; = 1.495 obtainecP^ simply by 
setting e = d u (l) — d — 3/2 in the two- loop e-expansion. 
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Starting with the HT expansions of £j_(Ki, R), similar considerations yield the estimate 
u± = 0.805 ± 0.015 ± 0.1\Rlp + 0.2733|. For this exponent, no MC results are available and 
our result can be compared only with the value u± = 0.757 from the e-expansion. 

The direct estimate of un from the analysis of £y is notoriously difficult because, in the 
P phase near the Lifshitz point, the competition between the nn and the nnn interaction 
reduces drastically the correlation length in the z direction. In particular, in the N — 1 
case, the determination of un has so far eluded even the most extensive^ MC simulation 
so far available. In the large N limit (namely in the case of the uniaxial spherical model), 
in which very long HT expansions can be easily computed, it was observed 25 that at least 
35 orders are necessary to approximate the behavior of £||! Also in the N = 2 case under 
study, the length of our HT expansions, unfortunately, seems to be still insufficient. We can, 
however, try to estimate indirectly un, either from a measure of the exponents f3 q and <f), 
taking advantage of the scaling law un = j3 q 4>, or by determining the exponent ai and then 
using the generalized hyperscaling law along with our previous estimate of u±. As above 
remarked, the determination of the exponent on from our series is not yet accurate enough 
to be useful. It does suggest, however, that a.\ = —0.02(2), which is compatible with the e- 
expansion estimate in Table [I] If, quite conservatively, we simply assume that a\ is negative 
and | ai| < 0.1, we get from the hyperscaling law the rough bounds 0.36 < u\\ < 0.52, which 
are consistent with the 0(e 2 ) value un ~ 0.372 reported in table [B On the other hand, using 
the above reported estimate = 1.00(4) and the estimate (3 q = 0.40(2) given below, we 
get the value un = 0.40(3). Plugging back this result into the hyperscaling law, we find the 
value ai = —0.01(6) for the specific-heat exponent, which is also compatible with all above 
indicated estimates. The values of the remaining exponents can all be obtained from the 
scaling laws. 

We can thus conclude that our HT determinations of the critical exponents at the LP 
are in most cases consistent, to within ~ 10%, with the e— expansion estimates. 

In order to map out the branch of the critical line with R < Rlp, we have to analyze the 
wave-vector-dependent susceptibility x(Qz, K\, R)- For R < Rlp the critical point Ki c (R) 
can be determined by locating the nearest positive singularity in K\ of x{Qz, Ki, R) with q z 
near the peak value q z = q z (Ki, R). The value of Ki c (R) is taken to be the minimum with 
respect to q z of the singularity locus Ki = Ki(q z , R) of x(Qz, K^R). Our results for the P-M 
branch of the critical line, also drawn in Fig{TJ complete the map of the boundary of the 
paramagnetic phase. In the same figure, we have schematically indicated a transition line, 
which separates the ferromagnetic and the modulated phases and joins the LP to the point 
R = — 1/4, T = 0, where the ordering of the ground state changes from ferromagnetic to 
modulated. This line is also expected to be of second order, but is beyond the reach of our 
HT analysis. Our phase diagram agrees well with the cited results^ obtained supplementing 
a MC simulation with the sixth-order HT expansions of Ref . (25| . 

It must be noted that the uncertainties of the points of the P-M branch of the critical 
line are sizably larger than those of the P-F branch, making it more difficult to obtain 
precise temperature-biased estimates of the critical exponents. If, nevertheless, we insist 
in computing some rough estimate of the exponent j(R), it is encouraging to observe that 
the results obtained from temperature-biased DAs remain essentially consistent, over a wide 
range of values of R, with those computed by the critical-point renormalization method 
which is insensitive to the uncertainties of the critical temperatures. The results for "f(R) 
obtained by these two methods are reported in Figj9j For R < —0.8, on the left of the region 
of crossover from the value of ji, we observe that both sequences of estimates tend to stabilize 
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at some value which is intermediate between those of the 0(2) and of the 0(3) universality 
classes. These puzzling results might deserve further confirmation by a MC study. If they 
are confirmed, the MC approach^! would also be best suited to investigate whether a possible 
(weak) first-order— (rather than second-order) character of the P-M transition might explain 
its features. This possibility was suggested by a renormalization group study of a Landau- 
Ginzburg effective four-component model of a biaxial N = 2 system- which turned out to 
have only unstable fixed points at second order in the e-expansion. 

The curve describing the peak value q z of the modulation wave-number q z at T = T C (R) 
vs R, obtained from our analysis, is reported in FigJTU] and compared to the mean-field 
prediction q h z IF {R) = cos' 1 (Rfjp / R) with i?f P F = -1/4. As pointed out in Ref.Q, also 
in the Ising case, at high temperature the peak of x(lz, K\, R) occurs at q^ IF (R) but, as 
the temperature decreases to T C (R), the peak moves to lower values of q z for R > —0.6 or 
otherwise to slightly higher values. This is clearly shown in FigJTUl The largest deviations 
of our results from the mean- value approximation occur in a small vicinity of Rlp- 

From the behavior of this curve near R ~ Rlp, we have estimated the value f3 q = 0.40(2) 
used above. 
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APPENDIX A: DERIVATION AND VALIDATION OF THE SERIES 

We have used the algorithm^ of the Schwinger- Dyson equations(SDE) to compute the 
HT expansions for the iV-vector spin models under study. When it was introduced, the 
SDE method could be profitably applied only to systems with nn interactions^, due to the 
severe limitations in the memory and speed of the computers available two decades ago. 
This method was repeatedly described in Refs. [70l , 71 1. and therefore it is sufficient to recall 



only that for the XY models described by the Hamiltonians ([I]) and (j2J), the SDE take a 
particularly simple and suggestive form. Indicating by s = (xi,qi,X2,q2',---Xk,qk) a set of 
site coordinates X{ and of integer nonzero quantities attached to them and such that 
Si=i Qi = 0' ^he generic correlation function can be written as 

< 0(s) >=\ J n^MCs)^* (Al) 

where 0(s) = exp[i J2k Qk&x k ] an d @x is the angle formed by the unit vector v{x) with a fixed 
direction. The SDE are: 

< m >= -— £(< 0«) > - < >) - — £(< 0(s+) > - < ^(s;) >) (A2) 
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Here 

00^) = exp(W*i) exp(-^x 1+ a M )0(s), (A3) 

while 

4>(%) = exp(-i0 Xl ) exp(i0 Xl+OM )0(s). (A4) 

Moreover 

</>(s+) = exp(i9 xi )exp(-i9 xi+bu )<f>(s), (A5) 

and 

0(s~) = exp(-^ Xl ) exp(ie xi+6 j0(s) (A6) 

with xi an arbitrary site in the set s. In the first sum of eq. flA2j) . a M is the lattice vector 
joining the site x\ to the nn sites, whereas in the second sum b v denotes the vector joining 
the site x\ to the nnn sites in the case of the Hamiltonian eq. (the vector joining the 
site X\ to the sn sites in the case of the Hamiltonian eq. (j2J)). 

The computational complexity of the SDE method increases with the lattice dimension- 
ality and, for a given dimensionality, with the effective coordination-number, i.e. as the 
number of interacting neighbors increases. Essentially the same difficulty is met with the 
evaluation of the embedding factors in the conventional graph approach, but our method 
has the advantage of using the SDE as recurrence relations among the correlation functions 
and of requiring only a straightforward iteration of these relations, which avoids altogether 
all combinatorial intricacies of the graph method. In the specific cases under study, due to 
the presence of the nnn interactions, this method is, by some orders of magnitude, more 
memory and CPU demanding than for the pure nn interactions, but otherwise not much 
more difficult. Thus for example, in the computationally most intensive case: that of the 
Hamiltonian of eq.([T]) in 3d with 3 — axial interaction, the general SDE is a linear relation- 
ship among 25 a priori different, correlation functions, while in the simpler nn interaction 
case in three- dimensions, the SDE involve only 13 correlation functions. 

A carefully designed code for the Schwinger-Dyson equations can compute moderately 
long series using only a reasonable CPU time of an ordinary 32-bit single-core processor 
desktop PC. In particular, our codes can reproduce all previously existing series data in 
a negligible CPU time of the order of 10~ 3 sec. Due to the sensitive dependence of the 
computational load on the effective coordination number, the series that we could derive for 
the 1 — axial models are longer (and more fastly computed) than those for the d — axial 
models. In the case of the uniaxial model studied here, the expansion through order 17 was 
completed in approximately three weeks by a PC, while the order 18 was obtained by using a 
few nodes of a PC cluster for an equivalent single-processor CPU-time of approximately five 
months. Longer series might be obtained by a more extensive parallelization of our codes 
or, perhaps, by returning to the conventional graph methods, provided that the efficiency 
of the current graph-embedding algorithms can be drastically improved. 

The comparison of the extended expansions with independent, sufficiently long and re- 
liable previous results is always a necessary step in the validation of the codes and of the 
results for automated derivations. Generally, in our case, this was not satisfactorily feasible. 
In two dimensions, a comparison was possible only for the 2-axial model (pQ), for which a 



fifth-order expansion of the susceptibility was tabulated in Ref . [40| and in three-dimensions, 
only in the case of the 1 -axial model, for which an expansion of the susceptibility was tabu- 
lated through sixth order in Ref.j25|. Our series agree with these results. Due to the lack of 
other published data and to the present unavailability of old short unpublished series^ for 
the 3-axial case in three dimensions, no further comparison with independent calculations 
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was possible. Weaker (and obvious) tests, through all orders that we have computed, are 
feasible in the limiting cases^ 3 - in which one of the exchange interaction constants vanishes. 
For example, when Ki = 0, the series for the m-axial models should reduce to those of the 
nn interaction models on the same lattice. Therefore the series for the 1- axial models should 
reduce to those of a nn interaction model on the same lattice when K 2 = 0, and to those of 
a nn model on a Id lattice when K 1 = 0. Of course, our results pass also these tests, which 
however pin down only 2 out of the I + 1 coefficients of each order /. Our confidence in the 
correctness of the calculations, for which we have written two sets of largely independent 
codes, in Fortran and in C++, receives further support also from the stability of the numer- 
ical results under our many rewritings of both sets of codes to increase their efficiency, as 
well as from the smooth and consistent behavior of the quantities analyzed. 
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FIG. 1: The phase diagram of the three-dimensional uniaxial XY model. In the (R,T C (R)) plane, 
we have represented by a continuous line the locus of critical points separating the disordered 
(paramagnetic) phase from the ordered (ferromagnetic and modulated) phases. In the scale of 
the figure, the uncertainties of the points are smaller than the width of the line. The Lifshitz 
point (LP) is located at the intersection of the critical and the disorder line (represented by a 
dashed curve). A transition line, also expected to be of second order, separates the modulated 
from the ferromagnetic phase, and joins the Lifshitz point to the point (R = — 1/4, T = 0), where 
the ordering of the ground state changes. This line cannot be mapped out by HT methods and 
therefore is only schematically indicated by a sequence of open circles. 
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FIG. 2: A representation of the modified-ratio sequences {K\ c (R)) n obtained from eg. H17H . For 
a few values of R chosen in a vicinity of Rm, the sequences have been extrapolated to large 
order by fitting them to the asymptotic form bi(R) — b2(R)/n 1+UJ with uj = 0.52. We have then 
plotted the modified-ratio sequences vs x = l/n 1+UJ , after normalizing them to their limiting 
values b\(R). The elements of the normalized sequences are represented by open triangles and 
are connected by continuous lines to guide the eye. The corresponding best fits to the asymptotic 
form: 1. — a(R)/n 1+ul with a{R) = 62 /b\{R) are represented by straight continuous lines. 



TABLE I: Values of the critical exponents for the (m, d, N) = (1, 3, 2) (uniaxial XY) Lifshitz point 
in three dimensions. The column with heading MF, taken from Ref. 50] reports the mean field 
exponents. In the column with heading 0(e 2 ), we report the values obtained^, by evaluating at 
e = 3/2 the e— expansions of the exponents truncated at the second order. Under the heading HT 
we report the results directly obtained in this paper. The heading MC refers to the results of the 
simulations in Refs. 4611471 ] . 



Exponent 


MF 


G(e 2 ) 


HT 


MC 




1 


1.495 


1.535(25) 


1.5(1) 


V\\ 


1 
4 


0.372 


0.40(3) 




Vl_ 


1 
2 


0.757 


0.805(15) 




V\\ 





-0.020 






V± 





0.042 






ai 





-0.047 




0.10(14)* 


A 


1 

2 


0.276 




0.20(2) 


Si 


3 








4>i 


1 
2 


0.725 


1.00(4) 






1 
2 


0.521 


0.40(2) 
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FIG. 3: The exponents "f(R) and u±(R) (notice that for R > u± = v) are computed along the 
P-F branch of the critical line K\ c = K\ C {R) in order to display their universality with respect to 
R. The horizontal continuous lines are bands of 0.5 % deviation from the central values, indicated 
by dashed lines, of recent high-precision estimates^ of the exponents 7 (in the case of the upper 
band), v (central band, shifted upwards by 0.6) and of the ratio v/j (lower band shifted upwards by 
0.7) for the XY universality class. We have indicated by open circles our estimates of the exponent 
j(R) obtained from DAs biased with the critical temperature, while the estimates obtained by the 
critical-point renormalization method, generally subject to a larger uncertainty, are denoted by 
open triangles. In the case of the exponent v^{R) we have denoted by open rhombs the estimates 
obtained from DAs biased with the critical temperature and by open squares the estimates obtained 
by the critical-point renormalization method. As already noticed, for graphical convenience, the 
values of this exponent are shifted upwards by 0.6. Finally a sequence of stars denotes the ratios 
u ±(R)/l(R) which, again for convenience, are shifted upwards by 0.7. A sequence of crosses 
schematically indicates a quantity proportional to the the absolute value |a+(i?)| of the correction- 
to-scaling amplitude. For graphical convenience, this quantity is shifted upwards by 1.4. The 
vertical dashed line on the left-hand side indicates the value Rlp = —0.2733(5) corresponding to 
the LP. 
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FIG. 4: A slightly modified and blown up representation of some of the data appearing in Figj3j 
We report the exponents 7(i?)(open circles), v±(R)(open squares) and the ratio u±(R) / j(R)(open 
triangles) computed by temperature-biased DAs along the P-F branch of the critical line K\ c = 
Ki c (R). The data are now normalized to the central values of the corresponding comparison 
estimates- 8 of the exponents for the XY universality class. The horizontal solid lines are bands of 
0.5 % deviation from the comparison value (central dashed line). The vertical dashed line on the 
left-hand side indicates the value of Rlp- The upper curve, schematically denoted by crosses, plots 
a quantity proportional to the the absolute value \a+(R)\ of the correction-to-scaling amplitude. 
This quantity is shifted upwards by 1.01 for graphical convenience. 



TABLE II: The critical values of K\ for selected values of R. The uncertainties reported here 
correspond only to the spread of the DA estimates and therefore are likely to underestimate the 
real errors when \R — Rm\ is not small, particularly so for R < 0. 



R 


K lc (R) 


R 


K lc (R) 


1.200 


0.15239(4) 


0.500 


0.18429(1) 


1.100 


0.15596(4) 


0.400 


0.19075(1) 


1.00 


0.15980(4) 


0.300 


0.19799(1) 


0.900 


0.16391(3) 


0.200 


0.20623(1) 


0.800 


0.16836(3) 


0.100 


0.21577(1) 


0.700 


0.17319(2) 


0.000 


0.22710(2) 


0.600 


0.17847(2) 


-0.100 


0.24113(4) 
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FIG. 5: Same as FigJH but here the exponents j(R) (open circles) and z^_i_(i?) (open squares) 
are computed along the P-F branch of the critical line by DAs biased with K\ C {R) + 5K\ C {R) /2 
for R> R M + 0.03 or biased with K Xc (R) - SK lc (R)/2 for R < R M - 0.03. Like in FiggJ the 
exponents are normalized to the mentioned 38 comparison central values 7 = 1.3178 and v = 0.6717, 
respectively. The horizontal solid lines are bands of 0.5 % deviation from the central value indicated 
by a dashed line. The vertical dashed line indicates the value of Rlp- The sequence of crosses 
represents a quantity proportional to the the absolute value |a^r(ii)| of the correction-to-scaling 
amplitude (shifted upwards by 1.01 for graphical convenience). 
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FIG. 6: Universality with respect to R of the moment ratios Q(p,q;r,s;R). We have reported 
results for Q(l/2, 1/2; 0, 1; R) represented by circles, for Q(l/2, 1/4; 0, 3/4; R) (triangles) and for 
Q(l/2, 1/2; 1/4, 3/4; i?) (rhombs). The ratios are normalized to their values at i?M- The values 
of these ratios at R M are respectively: Q(l/2, 1/2; 0, 1; R M ) = 0.876(1), Q(l/2, 1/4; 0, 3/4; i? M ) = 
0.929(1) and Q(l/2, 1/2; 1/4, 3/4; i? m ) = 0.9695(1). The sequence of crosses represents a quantity 
proportional to the the absolute value |a+(i?)| of the correction-to-scaling amplitude (and is shifted 
upwards by 1.01 for graphical convenience). The vertical dashed line indicates the value of Rlp- 
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FIG. 7: The correlation function C(0, 0, 0; 0, 0, d) between the spin at the origin and a spin on 
the function of the distance d between the spins, at fixed values of R and T. We have 

chosen the value {R = 0.11, T c (R)/2 = 2.59) (open circles) on the right-hand of the disorder line 
and the value (R = — 1.0, T c (R)/2 = 2.439)(open triangles) on the left-hand of it, to display the 
different behavior of the corresponding correlation function (pure exponential decay and oscillating 
exponential decay respectively). The lines connecting the symbols are drawn as a guide for the 
eye. 
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FIG. 8: The energy (upper curve), the ron-spin correlation function (middle curve) and the nnn- 
spin correlation function (lowest curve) computed near the boundary of the paramagnetic phase 
(at T = 1.1T C (R)) and plotted vs R. The vertical dashed line indicates the value of Rlp- 
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FIG. 9: The exponent 7(-R) computed along the P-M branch of the critical line. The estimates 
obtained by the critical-point renormalization method are represented as open triangles, those 
obtained by DAs biased with the critical temperature as open circles. The horizontal continuous 
lines represent the values of the exponent 7 for the 0(2) (lowest line), 0(3) (intermediate line) and 
0(4) (upper line) universality classes. The vertical dashed line indicates the value of Rlp- 
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FIG. 10: The peak value of the modulation vector q z at the critical temperature plotted vs. R for 
R < Rlp- The continuous curve represents the mean-field prediction q^ F = cos^ 1 (R^p / R) with 
R%f = -1/4. 
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